Skewed Cauchy distribution (skewcauchy)

Skewed Cauchy distribution (skewcauchy)#

The skewcauchy distribution is a heavy-tailed continuous distribution that generalizes the Cauchy by introducing a skewness parameter a ∈ (-1, 1).

  • When a = 0, it reduces to the standard (symmetric) Cauchy.

  • For a > 0, the distribution has more probability mass on the right and a heavier right tail (large positive outliers are more common).

  • For a < 0, the left tail is heavier.

A key feature (shared with the Cauchy): the mean and variance do not exist.

Learning goals#

  • Understand the piecewise nature of the PDF/CDF and how a controls asymmetry.

  • Know which quantities are well-defined (quantiles, entropy) and which are not (mean/variance).

  • Implement NumPy-only sampling via the inverse CDF (and an equivalent mixture view).

  • Visualize PDF/CDF and verify Monte Carlo simulations.

  • Use scipy.stats.skewcauchy for evaluation, simulation, and MLE fitting.

import platform

import numpy as np

import plotly.graph_objects as go
import os
import plotly.io as pio
from plotly.subplots import make_subplots

import scipy
from scipy import special
from scipy.stats import binomtest, skewcauchy

import plotly

pio.templates.default = "plotly_white"
pio.renderers.default = os.environ.get("PLOTLY_RENDERER", "notebook")

np.set_printoptions(precision=4, suppress=True)
rng = np.random.default_rng(7)

print("Python", platform.python_version())
print("NumPy", np.__version__)
print("SciPy", scipy.__version__)
print("Plotly", plotly.__version__)

Python 3.12.9
NumPy 1.26.2
SciPy 1.15.0
Plotly 6.5.2

1) Title & classification#

  • Name: skewcauchy

  • Type: continuous distribution

  • Support: \(x \in (-\infty, \infty)\)

  • Parameter space:

    • shape / skewness: \(a \in (-1,1)\)

    • location: \(\text{loc} \in \mathbb{R}\)

    • scale: \(\text{scale} > 0\)

We write the standardized form as:

\[Y \sim \mathrm{SkewCauchy}(a)\]

and the location/scale family as:

\[X = \text{loc} + \text{scale}\,Y.\]

SciPy uses this parameterization: skewcauchy(a, loc=..., scale=...).

2) Intuition & motivation#

What it models#

skewcauchy models real-valued data with extremely heavy tails (Cauchy-like) and asymmetry:

  • you expect outliers to be common;

  • you expect those outliers to be more common on one side.

This is useful when Gaussian or even Student-\(t\) noise is too light-tailed, and when symmetry is not appropriate.

Typical use cases#

  • Robust residuals with skew: regression or measurement error where large deviations happen, but predominantly positive or predominantly negative.

  • Asymmetric shock models: systems where positive jumps are more likely than negative jumps (or vice versa).

  • A building block in mixture models when you want components that can be both heavy-tailed and skewed.

Relations to other distributions#

  • Cauchy: skewcauchy(a=0) is the usual Cauchy.

  • Half-Cauchy: conditional on the sign, the magnitude is half-Cauchy with a scale that depends on a.

  • Skewed generalized t: skewcauchy is a special case (degrees of freedom \(\nu=1\)) of a broader skewed-\(t\) family.

A convenient generative picture (standardized case):

  1. Draw a sign \(S\in\{-1,+1\}\) with \(\mathbb{P}(S=+1)=(1+a)/2\).

  2. Draw a magnitude \(R\ge 0\) from a half-Cauchy whose scale is \((1+a)\) if \(S=+1\) and \((1-a)\) if \(S=-1\).

  3. Return \(Y = S\,R\).

This makes the meaning of a very concrete: it controls both which side is more likely and how heavy the tail is on that side.

3) Formal definition#

Throughout this section, \(-1<a<1\).

PDF#

The standardized PDF used by SciPy is:

\[ f(y; a) = \frac{1}{\pi\left(\frac{y^2}{\left(1+a\,\mathrm{sign}(y)\right)^2} + 1\right)}. \]

Equivalently, it is piecewise:

\[\begin{split} f(y;a) = \begin{cases} \dfrac{1}{\pi\left(1 + \left(\dfrac{y}{1-a}\right)^2\right)}, & y<0,\\ \dfrac{1}{\pi\left(1 + \left(\dfrac{y}{1+a}\right)^2\right)}, & y\ge 0. \end{cases} \end{split}\]

The location/scale form is:

\[ f_X(x; a,\text{loc},\text{scale}) = \frac{1}{\text{scale}}\, f\!\left(\frac{x-\text{loc}}{\text{scale}}; a\right). \]

CDF#

The standardized CDF has a closed form:

\[\begin{split} F(y;a) = \begin{cases} \dfrac{1-a}{2} + \dfrac{1-a}{\pi}\arctan\left(\dfrac{y}{1-a}\right), & y\le 0,\\ \dfrac{1-a}{2} + \dfrac{1+a}{\pi}\arctan\left(\dfrac{y}{1+a}\right), & y>0. \end{cases} \end{split}\]

And for \(X = \text{loc}+\text{scale}\,Y\):

\[F_X(x)=F\!\left(\frac{x-\text{loc}}{\text{scale}};a\right).\]

We’ll implement pdf, cdf, ppf (inverse CDF), and sampling from scratch using only NumPy.

def _check_skewcauchy_params(a: float, scale: float) -> None:
    if not (-1.0 < float(a) < 1.0):
        raise ValueError("shape parameter a must satisfy -1 < a < 1")
    if float(scale) <= 0.0:
        raise ValueError("scale must be > 0")


def skewcauchy_pdf(x: np.ndarray, a: float, loc: float = 0.0, scale: float = 1.0) -> np.ndarray:
    _check_skewcauchy_params(a, scale)
    x = np.asarray(x, dtype=float)
    y = (x - loc) / scale

    s = np.sign(y)
    denom = np.pi * (y * y / (1.0 + a * s) ** 2 + 1.0)
    return 1.0 / (scale * denom)


def skewcauchy_logpdf(x: np.ndarray, a: float, loc: float = 0.0, scale: float = 1.0) -> np.ndarray:
    _check_skewcauchy_params(a, scale)
    x = np.asarray(x, dtype=float)
    y = (x - loc) / scale

    s = np.sign(y)
    return -np.log(scale) - np.log(np.pi) - np.log(y * y / (1.0 + a * s) ** 2 + 1.0)


def skewcauchy_cdf(x: np.ndarray, a: float, loc: float = 0.0, scale: float = 1.0) -> np.ndarray:
    _check_skewcauchy_params(a, scale)
    x = np.asarray(x, dtype=float)
    y = (x - loc) / scale

    left = (1.0 - a) / 2.0 + (1.0 - a) / np.pi * np.arctan(y / (1.0 - a))
    right = (1.0 - a) / 2.0 + (1.0 + a) / np.pi * np.arctan(y / (1.0 + a))
    return np.where(y <= 0.0, left, right)


def skewcauchy_ppf(q: np.ndarray, a: float, loc: float = 0.0, scale: float = 1.0) -> np.ndarray:
    _check_skewcauchy_params(a, scale)
    q = np.asarray(q, dtype=float)

    q0 = (1.0 - a) / 2.0  # F(0)
    left = (1.0 - a) * np.tan(np.pi / (1.0 - a) * (q - q0))
    right = (1.0 + a) * np.tan(np.pi / (1.0 + a) * (q - q0))

    y = np.where(q < q0, left, right)
    return loc + scale * y


def skewcauchy_rvs(
    a: float,
    loc: float = 0.0,
    scale: float = 1.0,
    size=None,
    *,
    rng=None,
    eps: float = 1e-12,
) -> np.ndarray:
    _check_skewcauchy_params(a, scale)
    rng = np.random.default_rng() if rng is None else rng

    u = rng.random(size)
    u = np.clip(u, eps, 1.0 - eps)
    return skewcauchy_ppf(u, a, loc=loc, scale=scale)


def skewcauchy_rvs_mixture(
    a: float,
    loc: float = 0.0,
    scale: float = 1.0,
    size=None,
    *,
    rng=None,
    eps: float = 1e-12,
) -> np.ndarray:
    _check_skewcauchy_params(a, scale)
    rng = np.random.default_rng() if rng is None else rng

    p_right = (1.0 + a) / 2.0
    u_sign = rng.random(size)
    s = np.where(u_sign < p_right, 1.0, -1.0)

    u = rng.random(size)
    u = np.clip(u, eps, 1.0 - eps)

    # Half-Cauchy inverse CDF: R = b * tan(pi * U / 2)
    b = np.where(s > 0, 1.0 + a, 1.0 - a)
    r = b * np.tan(0.5 * np.pi * u)

    y = s * r
    return loc + scale * y

# Quick numerical checks against SciPy

a = 0.4
loc, scale = -1.25, 2.0

x = rng.normal(size=15) * 3

pdf_err = np.max(np.abs(skewcauchy_pdf(x, a, loc, scale) - skewcauchy.pdf(x, a, loc=loc, scale=scale)))
cdf_err = np.max(np.abs(skewcauchy_cdf(x, a, loc, scale) - skewcauchy.cdf(x, a, loc=loc, scale=scale)))

qs = np.array([1e-6, 0.1, 0.5, 0.9, 1 - 1e-6])
ppf_err = np.max(np.abs(skewcauchy_ppf(qs, a, loc, scale) - skewcauchy.ppf(qs, a, loc=loc, scale=scale)))

print("max |pdf diff|:", pdf_err)
print("max |cdf diff|:", cdf_err)
print("max |ppf diff|:", ppf_err)

max |pdf diff|: 0.0
max |cdf diff|: 0.0
max |ppf diff|: 0.0

4) Moments & properties#

Mean, variance, skewness, kurtosis#

Like the Cauchy, skewcauchy is heavy-tailed enough that the usual moments do not exist:

  • Mean: does not exist (not finite)

  • Variance: does not exist (infinite)

  • Skewness: undefined

  • Kurtosis: undefined

SciPy reflects this by returning nan for stats(..., moments='mvsk').

MGF and characteristic function#

  • The MGF \(M(t)=\mathbb{E}[e^{tX}]\) does not exist for any \(t\ne 0\).

  • The characteristic function \(\varphi(t)=\mathbb{E}[e^{itX}]\) does exist for all real \(t\).

For the standardized \(Y\sim\mathrm{SkewCauchy}(a)\), the real part has a clean form:

\[\Re\,\varphi_Y(t) = \frac{1-a}{2}e^{-(1-a)|t|} + \frac{1+a}{2}e^{-(1+a)|t|}.\]

The imaginary part is nonzero when \(a\ne 0\) and can be written using special functions (hyperbolic sine/cosine integrals, Shi/Chi). We’ll compute it in code.

Entropy#

The differential entropy of the standardized distribution is

\[h(Y)=\log(4\pi),\]

and for \(X=\text{loc}+\text{scale}\,Y\):

\[h(X)=\log(4\pi\,\text{scale}).\]

Notably, it is independent of the skewness parameter a.

def skewcauchy_cf_standard(t: np.ndarray, a: float) -> np.ndarray:
    if not (-1.0 < float(a) < 1.0):
        raise ValueError("-1 < a < 1 required")

    t = np.asarray(t, dtype=float)
    u = np.abs(t)
    sgn = np.sign(t)

    b_minus = 1.0 - a
    b_plus = 1.0 + a

    re = 0.5 * b_minus * np.exp(-b_minus * u) + 0.5 * b_plus * np.exp(-b_plus * u)

    im = np.zeros_like(re)
    mask = u > 0

    if np.any(mask):
        um = u[mask]

        # I_s(u, b) = ∫_0^∞ sin(u x) / (x^2 + b^2) dx
        #           = (cosh(bu) Shi(bu) - sinh(bu) Chi(bu)) / b, for u>0, b>0.
        Shi_p, Chi_p = special.shichi(b_plus * um)
        Shi_m, Chi_m = special.shichi(b_minus * um)

        Isp = (np.cosh(b_plus * um) * Shi_p - np.sinh(b_plus * um) * Chi_p) / b_plus
        Ism = (np.cosh(b_minus * um) * Shi_m - np.sinh(b_minus * um) * Chi_m) / b_minus

        im[mask] = sgn[mask] * (b_plus**2 / np.pi * Isp - b_minus**2 / np.pi * Ism)

    return re + 1j * im


ts = np.linspace(-10, 10, 2001)
a_vals = [-0.7, 0.0, 0.7]

fig = make_subplots(rows=2, cols=1, shared_xaxes=True, subplot_titles=("Re φ(t)", "Im φ(t)"))

for a in a_vals:
    phi = skewcauchy_cf_standard(ts, a)
    fig.add_trace(go.Scatter(x=ts, y=np.real(phi), mode="lines", name=f"a={a:+.1f}"), row=1, col=1)
    fig.add_trace(go.Scatter(x=ts, y=np.imag(phi), mode="lines", name=f"a={a:+.1f}"), row=2, col=1)

fig.update_layout(height=600, title="Characteristic function of standardized skewcauchy")
fig.update_xaxes(title_text="t", row=2, col=1)
fig.show()

# Entropy is constant in a (standardized case)

for a in [-0.9, -0.5, 0.0, 0.5, 0.9]:
    print(f"a={a:+.1f} -> entropy {skewcauchy.entropy(a):.12f}")

print("log(4π) =", float(np.log(4 * np.pi)))

a=-0.9 -> entropy 2.531024246969
a=-0.5 -> entropy 2.531024246969
a=+0.0 -> entropy 2.531024246969
a=+0.5 -> entropy 2.531024246969
a=+0.9 -> entropy 2.531024246969
log(4π) = 2.5310242469692907

5) Parameter interpretation#

Shape / skewness parameter a#

A simple, very interpretable identity for the standardized distribution is:

\[\mathbb{P}(Y>0)=\frac{1+a}{2},\qquad \mathbb{P}(Y<0)=\frac{1-a}{2}.\]

So a is literally the difference in sign probabilities:

\[a = \mathbb{P}(Y>0) - \mathbb{P}(Y<0).\]

It also controls the tail constants:

  • as \(y\to +\infty\), \(f(y)\sim \dfrac{(1+a)^2}{\pi y^2}\)

  • as \(y\to -\infty\), \(f(y)\sim \dfrac{(1-a)^2}{\pi y^2}\)

So a>0 means a heavier right tail, and a<0 a heavier left tail.

Location loc and scale scale#

  • loc shifts the distribution: \(X=\text{loc}+\text{scale}\,Y\).

  • scale stretches it, and increases entropy by \(\log(\text{scale})\).

Important nuance: for skewcauchy, loc is not the median unless a=0.

# How the median changes with a

a_grid = np.linspace(-0.95, 0.95, 401)
med = skewcauchy_ppf(0.5, a_grid)  # standardized median
p_right = (1 + a_grid) / 2

fig = make_subplots(rows=1, cols=2, subplot_titles=("Median of Y", "P(Y>0)"))
fig.add_trace(go.Scatter(x=a_grid, y=med, mode="lines", name="median"), row=1, col=1)
fig.add_trace(go.Scatter(x=a_grid, y=p_right, mode="lines", name="P(Y>0)"), row=1, col=2)
fig.update_xaxes(title_text="a", row=1, col=1)
fig.update_xaxes(title_text="a", row=1, col=2)
fig.update_yaxes(title_text="median", row=1, col=1)
fig.update_yaxes(title_text="probability", row=1, col=2)
fig.update_layout(height=350, title="Parameter interpretation")
fig.show()

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Cell In[6], line 4
      1 # How the median changes with a
      3 a_grid = np.linspace(-0.95, 0.95, 401)
----> 4 med = skewcauchy_ppf(0.5, a_grid)  # standardized median
      5 p_right = (1 + a_grid) / 2
      7 fig = make_subplots(rows=1, cols=2, subplot_titles=("Median of Y", "P(Y>0)"))

Cell In[2], line 38, in skewcauchy_ppf(q, a, loc, scale)
     37 def skewcauchy_ppf(q: np.ndarray, a: float, loc: float = 0.0, scale: float = 1.0) -> np.ndarray:
---> 38     _check_skewcauchy_params(a, scale)
     39     q = np.asarray(q, dtype=float)
     41     q0 = (1.0 - a) / 2.0  # F(0)

Cell In[2], line 2, in _check_skewcauchy_params(a, scale)
      1 def _check_skewcauchy_params(a: float, scale: float) -> None:
----> 2     if not (-1.0 < float(a) < 1.0):
      3         raise ValueError("shape parameter a must satisfy -1 < a < 1")
      4     if float(scale) <= 0.0:

TypeError: only length-1 arrays can be converted to Python scalars

6) Derivations#

We’ll use the standardized PDF and set:

\[b_+ = 1+a,\qquad b_- = 1-a.\]

For \(y\ge 0\):

\[f(y)=\frac{b_+^2}{\pi(y^2+b_+^2)},\]

and for \(y<0\):

\[f(y)=\frac{b_-^2}{\pi(y^2+b_-^2)}.\]

Expectation (why it does not exist)#

The mean exists (as a Lebesgue integral) iff \(\mathbb{E}[|Y|]<\infty\).

Consider the truncated absolute first moment:

\[\mathbb{E}[|Y|\,\mathbf{1}\{|Y|\le A\}] = \int_{-A}^{A} |y| f(y)\,dy.\]

Using the piecewise form and symmetry of \(|y|\):

\[ \int_0^A \frac{y\,b^2}{\pi(y^2+b^2)}\,dy = \frac{b^2}{2\pi}\log\left(\frac{A^2+b^2}{b^2}\right). \]

So

\[ \mathbb{E}[|Y|\,\mathbf{1}\{|Y|\le A\}] = \frac{b_+^2}{2\pi}\log\left(\frac{A^2+b_+^2}{b_+^2}\right) + \frac{b_-^2}{2\pi}\log\left(\frac{A^2+b_-^2}{b_-^2}\right) \xrightarrow[A\to\infty]{} \infty. \]

Therefore \(\mathbb{E}[|Y|]=\infty\) and the mean is undefined.

A subtle but important extra point: for \(a\ne 0\), even the principal value mean diverges because the positive and negative tails have different constants.

Variance (why it does not exist)#

The truncated second moment is

\[\mathbb{E}[Y^2\,\mathbf{1}\{|Y|\le A\}] = \int_{-A}^{A} y^2 f(y)\,dy.\]

For \(y\ge 0\):

\[ \int_0^A \frac{y^2 b^2}{\pi(y^2+b^2)}\,dy = \frac{b^2}{\pi}\left(A - b\arctan\left(\frac{A}{b}\right)\right), \]

which grows linearly in \(A\). Hence \(\mathbb{E}[Y^2]=\infty\) and the variance does not exist.

Likelihood#

Given observations \(x_1,\dots,x_n\), the log-likelihood is

\[ \ell(a,\text{loc},\text{scale}) = \sum_{i=1}^n \Bigl[-\log(\text{scale}) - \log\pi - \log\bigl(1 + y_i^2/(1+a\,\mathrm{sign}(y_i))^2\bigr)\Bigr] \]

where \(y_i=(x_i-\text{loc})/\text{scale}\).

We can maximize this numerically (SciPy does this in skewcauchy.fit).

def trunc_abs_moment(A: np.ndarray, a: float) -> np.ndarray:
    A = np.asarray(A, dtype=float)
    b_plus, b_minus = 1 + a, 1 - a
    term_p = (b_plus**2) / (2 * np.pi) * np.log((A * A + b_plus**2) / (b_plus**2))
    term_m = (b_minus**2) / (2 * np.pi) * np.log((A * A + b_minus**2) / (b_minus**2))
    return term_p + term_m


def trunc_pv_mean(A: np.ndarray, a: float) -> np.ndarray:
    A = np.asarray(A, dtype=float)
    b_plus, b_minus = 1 + a, 1 - a
    term_p = (b_plus**2) / (2 * np.pi) * np.log((A * A + b_plus**2) / (b_plus**2))
    term_m = (b_minus**2) / (2 * np.pi) * np.log((A * A + b_minus**2) / (b_minus**2))
    return term_p - term_m


def trunc_second_moment(A: np.ndarray, a: float) -> np.ndarray:
    A = np.asarray(A, dtype=float)
    b_plus, b_minus = 1 + a, 1 - a
    term_p = (b_plus**2) / np.pi * (A - b_plus * np.arctan(A / b_plus))
    term_m = (b_minus**2) / np.pi * (A - b_minus * np.arctan(A / b_minus))
    return term_p + term_m


A = np.logspace(-1, 3, 300)

fig = make_subplots(
    rows=1,
    cols=3,
    subplot_titles=("E[|Y| 1{|Y|≤A}]", "E[Y 1{|Y|≤A}]", "E[Y^2 1{|Y|≤A}]"),
)

for a in [-0.7, 0.0, 0.7]:
    fig.add_trace(go.Scatter(x=A, y=trunc_abs_moment(A, a), mode="lines", name=f"a={a:+.1f}"), row=1, col=1)
    fig.add_trace(go.Scatter(x=A, y=trunc_pv_mean(A, a), mode="lines", name=f"a={a:+.1f}"), row=1, col=2)
    fig.add_trace(go.Scatter(x=A, y=trunc_second_moment(A, a), mode="lines", name=f"a={a:+.1f}"), row=1, col=3)

fig.update_xaxes(type="log")
fig.update_layout(height=350, title="Divergence of truncated moments")
fig.show()

7) Sampling & simulation#

NumPy-only inverse CDF sampling#

Because we have a closed-form CDF and PPF, inverse transform sampling is straightforward:

  1. Draw \(U\sim\mathrm{Unif}(0,1)\).

  2. Return \(Y = F^{-1}(U)\).

We must clip \(U\) away from \(0\) and \(1\) because the tails are so heavy that the tangent in the PPF can overflow.

Equivalent mixture sampler#

Using the mixture intuition from Section 2:

  1. Draw \(S\in\{-1,+1\}\) with \(\mathbb{P}(S=+1)=(1+a)/2\).

  2. Draw \(R\sim \mathrm{HalfCauchy}(\text{scale}=1+a)\) if \(S=+1\) else \(R\sim \mathrm{HalfCauchy}(\text{scale}=1-a)\).

  3. Return \(Y=S\,R\).

We’ll implement both and check they agree.

a = 0.6

s1 = skewcauchy_rvs(a, size=200_000, rng=rng)
s2 = skewcauchy_rvs_mixture(a, size=200_000, rng=rng)

qs = [0.01, 0.05, 0.25, 0.5, 0.75, 0.95, 0.99]
print("Quantiles (inverse-CDF sampler):", np.quantile(s1, qs))
print("Quantiles (mixture sampler):    ", np.quantile(s2, qs))

print("Empirical P(Y>0):", float(np.mean(s1 > 0)), "theory:", (1 + a) / 2)

Quantiles (inverse-CDF sampler): [-4.9974 -0.9494  0.1595  1.0747  2.9982 15.986  77.1294]
Quantiles (mixture sampler):     [-5.0484 -0.9641  0.1591  1.0691  2.9816 16.3031 82.2879]
Empirical P(Y>0): 0.80144 theory: 0.8

8) Visualization#

We’ll visualize:

  • the PDF for several a

  • the CDF for several a

  • Monte Carlo samples (histogram + empirical CDF)

To keep plots readable, we focus on a finite \(x\)-range; the true distribution has extremely large outliers.

x = np.linspace(-10, 10, 3000)
a_vals = [-0.8, -0.4, 0.0, 0.4, 0.8]

fig = make_subplots(rows=1, cols=2, subplot_titles=("PDF", "CDF"))

for a in a_vals:
    fig.add_trace(go.Scatter(x=x, y=skewcauchy_pdf(x, a), mode="lines", name=f"a={a:+.1f}"), row=1, col=1)
    fig.add_trace(
        go.Scatter(x=x, y=skewcauchy_cdf(x, a), mode="lines", name=f"a={a:+.1f}", showlegend=False),
        row=1,
        col=2,
    )

fig.update_xaxes(title_text="x", row=1, col=1)
fig.update_xaxes(title_text="x", row=1, col=2)
fig.update_yaxes(title_text="density", row=1, col=1)
fig.update_yaxes(title_text="probability", row=1, col=2)
fig.update_layout(height=450, title="skewcauchy PDF/CDF (standardized)")
fig.show()

a = 0.6
samples = skewcauchy_rvs(a, size=80_000, rng=rng)

x_min, x_max = -15, 20

x_grid = np.linspace(x_min, x_max, 1200)
emp_cdf = np.searchsorted(np.sort(samples), x_grid, side="right") / samples.size
th_cdf = skewcauchy_cdf(x_grid, a)

fig = make_subplots(rows=1, cols=2, subplot_titles=("Histogram + PDF", "Empirical CDF vs CDF"))

fig.add_trace(
    go.Histogram(
        x=np.clip(samples, x_min, x_max),
        nbinsx=120,
        histnorm="probability density",
        name="samples (clipped for display)",
        opacity=0.35,
    ),
    row=1,
    col=1,
)
fig.add_trace(go.Scatter(x=x_grid, y=skewcauchy_pdf(x_grid, a), mode="lines", name="PDF"), row=1, col=1)

fig.add_trace(go.Scatter(x=x_grid, y=emp_cdf, mode="lines", name="empirical"), row=1, col=2)
fig.add_trace(go.Scatter(x=x_grid, y=th_cdf, mode="lines", name="theoretical"), row=1, col=2)

fig.update_xaxes(title_text="x", row=1, col=1)
fig.update_xaxes(title_text="x", row=1, col=2)
fig.update_yaxes(title_text="density", row=1, col=1)
fig.update_yaxes(title_text="probability", row=1, col=2)
fig.update_layout(height=450, title=f"Monte Carlo check (a={a})")
fig.show()

9) SciPy integration (scipy.stats.skewcauchy)#

SciPy exposes the distribution as scipy.stats.skewcauchy.

Key methods:

  • skewcauchy.pdf(x, a, loc, scale)

  • skewcauchy.cdf(x, a, loc, scale)

  • skewcauchy.rvs(a, loc, scale, size, random_state)

  • skewcauchy.fit(data) (MLE)

Note: skewcauchy.stats(a, moments='mvsk') returns nans because moments do not exist.

# Frozen distribution and MLE fit demo

a_true, loc_true, scale_true = 0.5, -1.0, 2.0
rv = skewcauchy(a_true, loc=loc_true, scale=scale_true)

data = rv.rvs(size=1500, random_state=7)

# Fit returns (a_hat, loc_hat, scale_hat)
a_hat, loc_hat, scale_hat = skewcauchy.fit(data)

print("true:", (a_true, loc_true, scale_true))
print("fit: ", (a_hat, loc_hat, scale_hat))

ll_true = float(np.sum(skewcauchy.logpdf(data, a_true, loc=loc_true, scale=scale_true)))
ll_fit = float(np.sum(skewcauchy.logpdf(data, a_hat, loc=loc_hat, scale=scale_hat)))
print("loglik(true) =", ll_true)
print("loglik(fit)  =", ll_fit)

print("stats(mvsk)  =", skewcauchy.stats(a_true, moments='mvsk'))

true: (0.5, -1.0, 2.0)
fit:  (0.4900946289203073, -1.0407488749013953, 2.019498790841202)
loglik(true) = -4840.023312243398
loglik(fit)  = -4839.383368018449
stats(mvsk)  = (nan, nan, nan, nan)

10) Statistical use cases#

Hypothesis testing#

Because the mean/variance do not exist, mean-based tests are not appropriate.

A simple skewness-related identity for the standardized distribution is:

\[\mathbb{P}(Y>0)=\frac{1+a}{2}.\]

So if the location is known (or you have a robust location estimate), you can test skewness via a binomial sign test:

  • \(H_0: a=0\) implies \(\mathbb{P}(Y>0)=1/2\).

This is not the most powerful test (it ignores magnitudes), but it is robust and fast.

Bayesian modeling#

  • As a likelihood for skewed heavy-tailed noise.

  • As a robust alternative to Gaussian noise when outliers are frequent and asymmetric.

Generative modeling#

  • As a heavy-tailed component in mixture models (e.g., clustering with outliers).

  • As a noise model in simulations where extreme events must be common and directionally biased.

# Example: sign test for a = 0 (assuming loc=0 is known)

n = 200

a_true = 0.5
samples = skewcauchy_rvs(a_true, size=n, rng=rng)  # standardized
k_pos = int(np.sum(samples > 0))

res = binomtest(k_pos, n=n, p=0.5, alternative="two-sided")

print("k positives:", k_pos, "out of", n)
print("a_hat_from_signs =", 2 * (k_pos / n) - 1)
print("p-value for H0: a=0 ->", res.pvalue)

k positives: 140 out of 200
a_hat_from_signs = 0.3999999999999999
p-value for H0: a=0 -> 1.507061557352626e-08

11) Pitfalls#

  • Invalid parameters: you must have -1 < a < 1 and scale > 0.

  • Moments do not exist: avoid mean/variance-based estimators, CLT intuition, and moment matching.

  • Tails are extreme: plots and histograms often need clipping or robust axis limits.

  • Numerical issues in sampling: the PPF uses tan(·), which blows up near 0 and 1; clip uniforms away from the endpoints.

  • Fitting can be tricky: MLE is sensitive to extreme points and can be flat in some directions; good initialization helps.

  • Use logpdf in products: multiplying many densities underflows; summing log-densities is stable.

12) Summary#

  • skewcauchy is a skewed, heavy-tailed generalization of the Cauchy with shape parameter a∈(-1,1).

  • It has no finite mean or variance (and thus no skewness/kurtosis moments).

  • The PDF/CDF are piecewise and have closed forms; sampling is easy via the inverse CDF.

  • The skewness parameter a controls both sign probability \(\mathbb{P}(Y>0)=(1+a)/2\) and tail heaviness.

  • SciPy provides a full implementation via scipy.stats.skewcauchy, including fit.

Contents